{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integration over Polytopes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "#### Extra dependencies : matplotlib (if using methods : plot_polytope and plot_polynomial) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sympy import sqrt\n",
    "from sympy.abc import x, y, z\n",
    "from sympy.geometry import *\n",
    "from sympy.integrals.intpoly import *\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Methods : "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### polytope_integrate(poly, expr, **kwargs)\n",
    "    Integrates polynomials over 2/3-Polytopes.\n",
    "\n",
    "    This function accepts the polytope in `poly` and the function in `expr` (uni/bi/trivariate polynomials are\n",
    "    implemented) and returns the exact integral of `expr` over `poly`.\n",
    "    \n",
    "    Parameters\n",
    "    ---------------------------------------\n",
    "    1. poly(Polygon) : 2/3-Polytope\n",
    "    2. expr(SymPy expression) : uni/bi-variate polynomial for 2-Polytope and uni/bi/tri-variate for 3-Polytope\n",
    "    \n",
    "    Optional Parameters\n",
    "    ---------------------------------------\n",
    "    1. clockwise(Boolean) : If user is not sure about orientation of vertices of the 2-Polytope and wants\n",
    "       to clockwise sort the points.\n",
    "    2. max_degree(Integer) : Maximum degree of any monomial of the input polynomial. This would require \n",
    "     \n",
    "   #### Examples :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "triangle = Polygon(Point(0,0), Point(1,1), Point(1,0))\n",
    "plot_polytope(triangle)\n",
    "print(\"Area of Triangle with vertices : (0,0), (1,1), (1,0) : \", polytope_integrate(triangle, 1))\n",
    "print(\"x*y integrated over Triangle with vertices : (0,0), (1,1), (1,0) : \", polytope_integrate(triangle, x*y),\"\\n\")\n",
    "\n",
    "hexagon = Polygon(Point(0, 0), Point(-sqrt(3) / 2, 0.5),\n",
    "                  Point(-sqrt(3) / 2, 3 / 2), Point(0, 2),\n",
    "                  Point(sqrt(3) / 2, 3 / 2), Point(sqrt(3) / 2, 0.5))\n",
    "plot_polytope(hexagon)\n",
    "print(\"Area of regular hexagon with unit side length  : \", polytope_integrate(hexagon, 1))\n",
    "print(\"x + y**2 integrated over regular hexagon with unit side length  : \", polytope_integrate(hexagon, x + y**2))\n",
    "\n",
    "polys = [1, x, y, x*y]\n",
    "print(\"1, x, y, x*y integrated over hexagon : \", polytope_integrate(hexagon, polys, max_degree=2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### main_integrate3d(expr, facets, vertices, hp_params)\n",
    "    Function to translate the problem of integrating uni/bi/tri-variate\n",
    "    polynomials over a 3-Polytope to integrating over its faces.\n",
    "    This is done using Generalized Stokes's Theorem and Euler's Theorem.\n",
    "    \n",
    "    Parameters\n",
    "    ------------------\n",
    "    1. expr : The input polynomial\n",
    "    2. facets : Faces of the 3-Polytope(expressed as indices of `vertices`)\n",
    "    3. vertices : Vertices that constitute the Polytope\n",
    "    4. hp_params : Hyperplane Parameters of the facets\n",
    "    \n",
    "   #### Examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cube = [[(0, 0, 0), (0, 0, 1), (0, 1, 0), (0, 1, 1), (1, 0, 0), (1, 0, 1), (1, 1, 0), (1, 1, 1)],\n",
    "         [2, 6, 7, 3], [3, 7, 5, 1], [7, 6, 4, 5], [1, 5, 4, 0], [3, 1, 0, 2], [0, 4, 6, 2]]\n",
    "vertices = cube[0]\n",
    "faces = cube[1:]\n",
    "hp_params = hyperplane_parameters(faces, vertices)\n",
    "main_integrate3d(1, faces, vertices, hp_params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### polygon_integrate(facet, index, facets, vertices, expr, degree)\n",
    "    Helper function to integrate the input uni/bi/trivariate polynomial\n",
    "    over a certain face of the 3-Polytope.\n",
    "    \n",
    "    Parameters\n",
    "    ------------------\n",
    "    facet : Particular face of the 3-Polytope over which `expr` is integrated\n",
    "    index : The index of `facet` in `facets`\n",
    "    facets : Faces of the 3-Polytope(expressed as indices of `vertices`)\n",
    "    vertices : Vertices that constitute the facet\n",
    "    expr : The input polynomial\n",
    "    degree : Degree of `expr`\n",
    "    \n",
    "   #### Examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cube = [[(0, 0, 0), (0, 0, 5), (0, 5, 0), (0, 5, 5), (5, 0, 0),\n",
    "                 (5, 0, 5), (5, 5, 0), (5, 5, 5)],\n",
    "                 [2, 6, 7, 3], [3, 7, 5, 1], [7, 6, 4, 5], [1, 5, 4, 0],\n",
    "                 [3, 1, 0, 2], [0, 4, 6, 2]]\n",
    "facet = cube[1]\n",
    "facets = cube[1:]\n",
    "vertices = cube[0]\n",
    "print(\"Area of polygon < [(0, 5, 0), (5, 5, 0), (5, 5, 5), (0, 5, 5)] > : \", polygon_integrate(facet, 0, facets, vertices, 1, 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### distance_to_side(point, line_seg)\n",
    "\n",
    "    Helper function to compute the distance between given 3D point and\n",
    "    a line segment.\n",
    "\n",
    "    Parameters\n",
    "    -----------------\n",
    "    point : 3D Point\n",
    "    line_seg : Line Segment\n",
    "    \n",
    "#### Examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "point = (0, 0, 0)\n",
    "distance_to_side(point, [(0, 0, 1), (0, 1, 0)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### lineseg_integrate(polygon, index, line_seg, expr, degree)\n",
    "    Helper function to compute the line integral of `expr` over `line_seg`\n",
    "\n",
    "    Parameters\n",
    "    -------------\n",
    "    polygon : Face of a 3-Polytope\n",
    "    index : index of line_seg in polygon\n",
    "    line_seg : Line Segment\n",
    "#### Examples :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "polygon = [(0, 5, 0), (5, 5, 0), (5, 5, 5), (0, 5, 5)]\n",
    "line_seg = [(0, 5, 0), (5, 5, 0)]\n",
    "print(lineseg_integrate(polygon, 0, line_seg, 1, 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### main_integrate(expr, facets, hp_params, max_degree=None)\n",
    "\n",
    "    Function to translate the problem of integrating univariate/bivariate\n",
    "    polynomials over a 2-Polytope to integrating over its boundary facets.\n",
    "    This is done using Generalized Stokes's Theorem and Euler's Theorem.\n",
    "\n",
    "    Parameters\n",
    "    --------------------\n",
    "    expr : The input polynomial\n",
    "    facets : Facets(Line Segments) of the 2-Polytope\n",
    "    hp_params : Hyperplane Parameters of the facets\n",
    "\n",
    "    Optional Parameters:\n",
    "    --------------------\n",
    "    max_degree : The maximum degree of any monomial of the input polynomial.\n",
    "    \n",
    "#### Examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "triangle = Polygon(Point(0, 3), Point(5, 3), Point(1, 1))\n",
    "facets = triangle.sides\n",
    "hp_params = hyperplane_parameters(triangle)\n",
    "print(main_integrate(x**2 + y**2, facets, hp_params))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### integration_reduction(facets, index, a, b, expr, dims, degree)\n",
    "    This is a helper function for polytope_integrate. It relates the result of the integral of a polynomial over a\n",
    "    d-dimensional entity to the result of the same integral of that polynomial over the (d - 1)-dimensional \n",
    "    facet[index].\n",
    "    \n",
    "    For the 2D case, surface integral --> line integrals --> evaluation of polynomial at vertices of line segments\n",
    "    For the 3D case, volume integral --> 2D use case\n",
    "    \n",
    "    The only minor limitation is that some lines of code are 2D specific, but that can be easily changed. Note that\n",
    "    this function is a helper one and works for a facet which bounds the polytope(i.e. the intersection point with the\n",
    "    other facets is required), not for an independent line.\n",
    "    \n",
    "    Parameters\n",
    "    ------------------\n",
    "    facets : List of facets that decide the region enclose by 2-Polytope.\n",
    "    index : The index of the facet with respect to which the integral is supposed to be found.\n",
    "    a, b : Hyperplane parameters corresponding to facets.\n",
    "    expr : Uni/Bi-variate Polynomial\n",
    "    dims : List of symbols denoting axes\n",
    "    degree : Degree of the homogeneous polynoimal(expr)\n",
    "    \n",
    "   #### Examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "facets = [Segment2D(Point(0, 0), Point(1, 1)), Segment2D(Point(1, 1), Point(1, 0)), Segment2D(Point(0, 0), Point(1, 0))]\n",
    "print(integration_reduction(facets, 0, (0, 1), 0, 1, [x, y], 0))\n",
    "print(integration_reduction(facets, 1, (0, 1), 0, 1, [x, y], 0))\n",
    "print(integration_reduction(facets, 2, (0, 1), 0, 1, [x, y], 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### hyperplane_parameters(poly) :\n",
    "    poly : 2-Polytope\n",
    "    \n",
    "    Returns the list of hyperplane parameters for facets of the polygon.\n",
    "    \n",
    "    Limitation : 2D specific.\n",
    "   #### Examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "triangle = Polygon(Point(0,0), Point(1,1), Point(1,0))\n",
    "hyperplane_parameters(triangle)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### best_origin(a, b, lineseg, expr) :\n",
    "    a, b : Line parameters of the line-segment\n",
    "    expr : Uni/Bi-variate polynomial\n",
    "    \n",
    "    Returns a point on the lineseg whose vector inner product with the divergence of expr yields an expression with \n",
    "    the least maximum total power. This is for reducing the number of computations in the integration reduction call.\n",
    "    \n",
    "    Limitation : 2D specific.\n",
    "    \n",
    "   #### Examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Best origin for x**3*y on x + y = 3 : \", best_origin((1,1), 3, Segment2D(Point(0, 3), Point(3, 0)), x**3*y))\n",
    "print(\"Best origin for x*y**3 on x + y = 3 : \",best_origin((1,1), 3, Segment2D(Point(0, 3), Point(3, 0)), x*y**3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### decompose(expr, separate=False) :\n",
    "    expr : Uni/Bi-variate polynomial.\n",
    "    separate(default : False) : If separate is True then return list of constituting monomials.\n",
    "    \n",
    "    Returns a dictionary of the terms having same total power. This is done to get homogeneous polynomials of\n",
    "    different degrees from the expression.\n",
    "    \n",
    "   #### Examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(decompose(1 + x + x**2 + x*y))\n",
    "print(decompose(x**2 + x + y + 1 + x**3 + x**2*y + y**4 + x**3*y + y**2*x**2))\n",
    "print(decompose(x**2 + x + y + 1 + x**3 + x**2*y + y**4 + x**3*y + y**2*x**2, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### norm(expr) :\n",
    "     \n",
    "    point : Tuple/SymPy Point object/Dictionary\n",
    "    \n",
    "    Returns Euclidean norm of the point object.\n",
    "\n",
    "   #### Examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(norm((1, 2)))\n",
    "print(norm(Point(1, 2)))\n",
    "print(norm({x: 3, y: 3, z: 1}))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### intersection(lineseg_1, lineseg_2) :\n",
    "     \n",
    "    lineseg_1, lineseg_2 : The input line segments whose intersection is to be found.\n",
    "    \n",
    "    Returns intersection point of two lines of which lineseg_1, lineseg_2 are part of. This function is\n",
    "    called for adjacent line segments so the intersection point is always present with line segment boundaries.\n",
    "\n",
    "   #### Examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(intersection(Segment2D(Point(0, 0), Point(2, 2)), Segment2D(Point(1, 0), Point(0, 1))))\n",
    "print(intersection(Segment2D(Point(2, 0), Point(2, 2)), Segment2D(Point(0, 0), Point(4, 4))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### is_vertex(ent) :\n",
    "     \n",
    "    ent : Geometrical entity to denote a vertex.\n",
    "    \n",
    "    Returns True if ent is a vertex. Currently tuples of length 2 or 3 and SymPy Point object are supported.\n",
    "   #### Examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(is_vertex(Point(2, 8)))\n",
    "print(is_vertex(Point(2, 8, 1)))\n",
    "print(is_vertex((1, 1)))\n",
    "print(is_vertex([2, 9]))\n",
    "print(is_vertex(Polygon(Point(0, 0), Point(1, 1), Point(1, 0))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### plot_polytope(poly) :\n",
    "     \n",
    "    poly : 2-Polytope\n",
    "    \n",
    "    Plots the 2-Polytope. Currently just defers it to plotting module in SymPy which in turn uses matplotlib.\n",
    "        \n",
    "   #### Examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hexagon = Polygon(Point(0, 0), Point(-sqrt(3) / 2, 0.5),\n",
    "                  Point(-sqrt(3) / 2, 3 / 2), Point(0, 2),\n",
    "                  Point(sqrt(3) / 2, 3 / 2), Point(sqrt(3) / 2, 0.5))\n",
    "plot_polytope(hexagon)\n",
    "\n",
    "twist = Polygon(Point(-1, 1), Point(0, 0), Point(1, 1), Point(1, -1),\n",
    "                Point(0, 0), Point(-1, -1))\n",
    "plot_polytope(twist)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### plot_polynomial(expr) :\n",
    "     \n",
    "    expr : The uni/bi-variate polynomial to plot\n",
    "    \n",
    "    Plots the polynomial. Currently just defers it to plotting module in SymPy which in turn uses matplotlib.\n",
    "        \n",
    "   #### Examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "expr = x**2\n",
    "plot_polynomial(expr)\n",
    "\n",
    "expr = x*y\n",
    "plot_polynomial(expr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
